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A theoretical analysis of Coulomb systems on lattices in general dimensions is 
presented. The thermodynamics is developed using Debye-Hiickel theory with ion- 
pairing and dipole-ion solvation, specific calculations being performed for 3D lattices. 
As for continuum electrolytes, low-density results for sc, bcc and fee lattices indicate 
the existence of gas- liquid phase separation. The predicted critical densities have 
values comparable to those of continuum ionic systems, while the critical tempera- 
tures are 60-70% higher. However, when the possibility of sublattice ordering as well 
as Debye screening is taken into account systematically, order-disorder transitions 
and a tricritical point are found on sc and bcc lattices, and gas-liquid coexistence is 
suppressed. Our results agree with recent Monte Carlo simulations of lattice elec- 
trolytes. 
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I. INTRODUCTION 

It is well known that criticality in simple fluids with short-range potentials can be de- 
scribed by the Ising universality class with critical exponents accessible via renormalization 
group calculations!. However, for Coulomb systems, where particles interact through long- 
range potentials, the nature of criticality remains open to question!. Numerous theoretical, 
experimental and computational investigations of electrolyte systems have not yet produced 
a clear picture of the thermodynamics in the critical region. Early experiments on criticality 
in electrolytes suggested a strong dichotomy: namely, some electrolytesi'S, termed solvopho- 
bic and typically having large solvent dielectric constant, are satisfactorily characterized by 
Ising critical exponents. This suggests, that the principal interactions driving the phase sep- 
aration in such systems are of short-range character. On the other hand, a number of organic 
salts in appropriate solvents, typically of low dielectric constant, were found to exhibit classi- 
cal or close-to-classical behavior!, and have been called Coulombic, stressing the importance 
of the dominant electrostatic interactions. Moreover, in sodium-ammonia solutions! (and 
some other systems: see!'! and references therein) , crossover from classical to Ising behavior 
had been observed, but at a reduced temperature t = (T — T c )/T c ~ 0.6 x 1CT 2 , unusually 
close to the critical point. This led to the ideal that the true asymptotic critical behavior 
of ionic fluids is always of Ising character but that crossover from nonasymptotic, close-to- 
classiea, behavior „ at S ca,e S t ba t ,aa y senses be exper— y inaeeessibl M. 

Monte Carlo computer simulations provide another useful method of investigating the 
properties of ionic systems. In the last decade substantial progress has been achieved in this 
field, with primary effort focused on the coexistence curvesi'00. However, some special 
attention has also been devoted to the heat capacity which is significant for elucidating the 
critical regioni00. Nevertheless, because of the long-range character of the interactions 
and the low values of the critical temperatures, which lead to many strongly bound ion 
pairsO'0, computer simulations for finite systems still lack the ability to clearly determine 
critical exponents and hence to identify the nature of the criticality. 

The success of the renormalization group (RG) approach in describing non-ionic fluidH 
suggests that it might also be applied to Coulombic criticality. However, to implement an 
RG treatment, the existence of a physically well based mean-field theory turns out to be 
crucial!. The simplest model for theoretical investigations of ionic systems is the so called 
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restricted primitive model (RPM), which considers particles of equal sizes and positive and 
negative charges of equal magnitude. Two main theoretical approaches have emerged. The 
first employs an extensioni'O'El'0 of the basic Debye-Hiickel (DH) theory!!, developed in 
the early 20th century for dilute solutions of strong electrolytes. The second approach rests 
on integral equations for correlation functions, typically employing the Ornstein-Zernike 
equations in combination with some truncation, as in the mean spherical approximation 
Neither of these two approaches has any known independent basis, such as 
an overall variational principle for the ionic free energy, that might help justify its reliability. 
However, compared to the values predicted by DH-based treatments, MSA-based theories 
yield relatively poor agreement with the critical parameters found by current Monte Carlo 
simulations, namely, in reduced unitsi'B, T* ~ 0.049 and p* = 0.06 to 0.085. Careful 
analysisiSi, utilizing thermodynamic energy bounds, etc., also suggests that DH-based 
theories promise a better description of the critical region of model electrolytes. 

Since the Ising model, which is equivalent to a lattice gas, has played a crucial role in un- 
derstanding critical phenomena in non-ionic systems, lattice models of electrolytes deserve 
attention. Although clearly artificial as regards the description of real ionic solutions, which 
possess continuous rather than discrete spatial symmetry, they are attractive for various 
reasons. First, by virtue of the lattice character one can incorporate the behavior of dense 
phases, at least one of which should be an ordered ionic crystal. Lattice models may also be 
effective for describing defects in real crystalsil. Second, even finely discretized lattice sys- 
tems present a computational advantage over their continuous-space counterparts in Monte 
Carlo simulations^. Last but not the least, discrete-state lattice models facilitate the deriva- 
tion of equivalent field-theoretical descriptions and, thereby, the study of the significance of 
various terms in the effective Hamiltonian. Moreover, Coulomb interactions can be exactly 
represented in terms of a nearest-neighbor Hamiltonian via the sine-Gordon transformation 
that yields both low fugacity and high-temperature expansions of the equation of state. (For 
a recent overview and results see Ref.0 and references therein.) 

Despite their distinct theoretical advantages, lattice models of electrolytes have not been 
studied systematically. The aim of the work reported below has been to repair this omission. 

Most of the previous analyticaSii anc i numerical0S@I!I3 work on lattice ionic systems 
has addressed the question of tricriticality and of order-disorder transitions. While the 
overall density is the order parameter which suffices to reveal gas-liquid critical behavior 
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in ionic solution, the presence of an underlying lattice allows naturally for the appearance 
of another order parameter. In bipartite lattices, such as the simple cubic (sc) and body- 
centered cubic (bcc) lattices, ions of opposite charges can distribute unequally between the 
sublattices, thereby reducing the electrostatic part of the free energy. At the same time, 
the entropy is also reduced which increase the free energy. This competition leads to the 
appearance of a phase with long-range order resembling an ionic crystal; second-order phase 
transitions are then a likely consequence. In the continuum case analogous oscillations 
appear in the charge-charge correlation functions^ at certain values of density-temperature 
ratio. However, such ordered phases may turn out to be thermodynamically so stable, that 
a gas-liquid phase transition predicted by a continuum theory may not survive in a lattice 
model: the lattice system tends to "solidify" before forming a "proper" liquid. Indeed, this 
scenario has been observed in numerical studies. On the other hand, the possible presence 
of both gas-liquid and tricritical points has been predicted theoretically by Ciach and Stell 
for a model with additional short-range interactions added to the lattice Coulomb forces0. 

As indicated, we present in this article a study of the simplest, single-site hard core model 
of lattice ionic system with charges q± = ±q. In Sec. II we describe the basic DH theory 
on general, (^-dimensional Bravais lattices. Our analysis focuses on d = 3 in Sec. III. After 
presenting the results for pure DH theory, the crucial phenomenon of Bjerrum ion pairing 
is introduced in Sec. III.B; but, following Fisher and LevinEl, this must be supplemented 
by explicit dipole-ion solvation effects: see Sec. III.C. Then, in Sec. IV the possibility of 
sublattice charge ordering is discussed. Unlike previous treatments^!!, we account for both 
electrostatic screening and sublattice ordering in a unified framework. Our conclusions are 
summarized briefly in Sec. V. 



II. LATTICE DEBYE-HUCKEL THEORY IN GENERAL DIMENSIONS 

Our derivation for general <i-dimensional lattices follows closely the Debye-Hiickel 
approachll. We confine ourselves to the lattice restricted primitive model (LRPM), which 
consists of oppositely charged ions with charges q and — q which occupy single lattice sites 
of a c/-dimensional Bravais lattice. In this simplest model the ions interact only through the 
electrostatic field and otherwise behave as ideal particles, subject only to on-site exclusion. 
Thus the total free energy density, which plays the central role in the thermodynamics of the 
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system, can be written / = f Id + f DH . As the overall system must be neutral, the average 
densities of the positive and negative ions are equal: p + = p_ = \p\. Correspondingly, for 
the reduced chemical potential and pressured we have 

p. + = p.-=p.i, P = max[/+//ipi], (1) 

pi 

where ft = p,/ksT and p = p/ksT. The ideal lattice gas contribution to the free energy is, 
up to a constant term, given by 

/" = ~ = ~ Pl M " X ~ ZA ln(l - Pi), (2) 

kTV Vq Vq 

with the corresponding chemical potential 

-u = _ d fid/ dpi = lnp * _ ln(1 _ p ^ (3) 

where V is the total lattice volume while pi = piv is the reduced dimensionless density of 
(free) ions and vq is the volume per site of the lattice. 

Next we determine the contribution to free energy arising from the Coulombic interac- 
tions. However, the lattice form of the potential, which takes into account the discreteness of 
the space, should be used. This lattice Coulomb potential will approach the continuous 1/r 
potential asymptotically at large distances, but it differs significantly at small distances. We 
start with the linearized lattice Poisson-Boltzmann equation, which determines the average 
electrostatic potential at point r. Following the standard DH approach^, we easily find 

A</?(r) = «VW - (qC d /Dv Q )5(r), (4) 

where k 2 = Cdfipiq 2 /D is the inverse Debye screening length, with (3 = 1/ksT. The constant 
factor Cd = 2n d / 2 /T(d/2) is determined by the dimensionality of the lattice system@. In 
this equation we use the lattice Laplacian defined through 

Ap(r) = ^ J2 ^ r + a ) - V^)} > (5) 
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where a is a nearest-neighbor vector and the summation runs over all Cq nearest neighbors. 
The DH equation, (|4]), can be easily solved by Fourier transformation yielding 

= Cd^a^ r 

V{r} D V{) 2d] k {xi + 2d)/2d-J{\Ly [) 



with x — Ka and f k = (2tt) d f_ d d \a. The lattice function J(k) is defined by 

J ( k ) = rE e *' a ( ? ) 

c o — 

nn 

In the DH approach, we need only the potential felt by an ion fixed at the origin due to all 
the surrounding ions. Because of the Bravais symmetry, we can find the total potential at 
the origin by averaging over the nearest-neighbor sites to obtain 



ip(r = 0) = — Vy?(a ran ). 



Co 

I bit 

Introducing the integrated lattice Green's function via 

1 



(8) 



and using (0) we obtain 



<P(0) 



p(0 



C d g a 2 
Dv 2d 



1 - C^(k) 



(9) 



P 



2d 



x 2 + 2d 



- 1 



(10) 



The potential due to a single ion in the absence of other ions is obtained simply by setting 
x — in (H) which yields 

C d q a 2 



<M0) 



(11) 



Dv 2d 

Subtracting this expression from the total potential (|10"D and using @ yields the potential 
felt by an ion of charge g« at r = due to all the surrounding ions, namely, 



C d qi a 2 
Dv 2d 



P 



2d 



P(D 



(12) 



x 2 + 2d, 

The electrostatic part of the free energy can be found by using the Debye charging 
procedure^! which yields 



El 



kTV 
1 



Advo 



1_ J2 F lon = -nfa 1 MM)dX 

d{x 2 



x 2 P(l) 



2d 



x 2 + 2d 



(13) 



Combining the ideal-gas and Debye-Hiickel terms yields / = / + / and 
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with reduced temperature denned by 

T* = kTDa d - 2 /q 2 . (15) 
From this we find the pressure for an arbitrary <i-dimensional Bravais lattice to be 

(16) 



pv = - ln(l + — 



Eqs.@ and (|T^ ) - ([16|) give full information about the thermodynamic behavior of a lattice 
Coulomb system. In particular, the possibility of phase transitions and criticality can be 
investigated by analyzing the spinodals, and the phase coexistence curves may be obtained 
by the matching pressure and chemical potential in coexisting phases. Spinodals are specified 
by setting the inverse isothermal compressibility to zero, so that 

PiiT- = °> ( 17 ) 



pikTKx dpi 



which, on taking (113) into account, reduces to 



C d a d gi-QdP(0/d( 
s 2dv 2 + (l-C) 2 dP(()/dC 1 J 

with ( = 2d/(x 2 + 2d). One can show that when p\ becomes large (which corresponds to 
C — > 0), one has T* « c Cda d (pl) 2 /2dv . Eqs.(|T^)-([T^) can be used to investigate the phase 
behavior of electrolytes in any dimension. 

We mention briefly here the critical parameters obtained for d = 1 and 2. In the one- 
dimensional case the lattice Green's function gives 



P(Q=ir dk = 1 as) 

ig Wo 1-CcosA; ^1^2' [U) 
which yields a spinodal of the form 



Ts * " x [x + (a* + 4)3/2] > x ~ Ka - ( 2 °) 
This specifies a critical point with parameters 

p* c = 0, T; = oo, (21) 

in accordance with the general principle that one-dimensional systems do not display phase 
transitions. However, since ip(r) oc |r| in a one-dimensional system, the DH method fails 
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Table I: Lattice parameters 



lattice 


unit cell 


nearest neighbor 


number of nearest 


volume per 




edge 


distance, a nn 


neighbors, cq 


site, Vo 


sc 


a 


a 


6 


"o 


bcc 


2a 


a V3 


8 




fee 


2a 


a V2 


12 


2a 3 



and describing the one-dimensional ionic lattice model demands a different approach: one 
may note the continuum analysis!!. 

For d = 2 dimensions the lattice Green's functions are given in RefM Then for both for 
triangular and square lattices the predicted critical parameters are found to be 

p* c = 0, T* c = 1/4, (22) 

precisely, the same values as for the continuum model^l. 



III. ELECTROLYTES IN THREE DIMENSIONS 



A. Pure DH theory 

Let us now examine 3D cubic lattices in more detail. We address three cases: simple 
cubic (sc), body centered cubic (bcc) and face centered cubic (fee); for convenience their 
geometrical parameters are listed in Table |. The basic lattice function J(k), defined in (^), 
is then given by 

J(k) = | (cos k\ + cos k 2 + cos kz), sc, (23) 

= cos/^i cos/c 2 cos/c 3 , bcc, (24) 

= | (cos ki cos k 2 + cos k 2 cos k 3 + cos k\ cos k 3 ), fee, (25) 

with — 7r < ki,k 2 ,k 3 < 7T. The corresponding integrated lattice Green's functions can be 
explicitly calculated using their representation in terms of complete elliptic integrals as 
shown by JoyceS The self-potential of an ion (in the absence of any screening) is then 
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given by 

(Da/q)MO) = (26) 
~ 1.082, 1.070 and 1.021, 

for sc, bcc and fee lattices. This reduced value approaches the exact continuum potential 
value 1 as the number of nearest neighbors increases. At low densities the free energy as 
given by flHf ) can be expanded in powers of x = na, which yields 

f El = TT^t 1 " 0.282X + 0.025x 2 + ...), sc, (27) 



x 



l- 0.286x + 0.025a; 2 + ...), bcc, (28) 



12vra 3 

;i- 0.296x + 0.025x 2 + ...), fee. (29) 



x 3 



12vra 3 

The leading term precisely reproduces the exact continuum DH result, which, of course, 
is independent of a. The magnitude of the first correction term increases with increasing 
coordination number; in the hard sphere continuum model it becomes 0.75. 

The predicted coexistence curves for the sc, bcc and fee lattices are shown in Fig. 1, while 
the critical parameters are listed in Table II. A surprising feature of these coexistence curves 
is that the liquid density approaches a finite value, p* iq (0), as T — ► that is substantially 
smaller than the maximum, close-packing density p\ = 1: see Table II. 

For comparison, Fig. 1 also displays the predictions of DH theory for the continuum 
RPM supplemented by hard-core interactions in the free volume approximation with the 
simple cubic packing limitS Although the critical temperatures decrease slightly as the 
number of nearest neighbors approaches realistic values of the coordination number, say, 
12 — 14, as observed in simple liquids, its value for all three lattices remains about 50% 
higher than the corresponding continuum value. This is, indeed, a rather general feature of 
lattice models, which tend to display higher critical temperatures than their continuum-space 
counterparts. However, the predicted critical densities are quite comparable, decreasing from 
about 60% above the continuum value to only 3 or 4 % higher: see Table III. Clearly, packing 
considerations play a significant role in the value of p*. 
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B. Bjerrum ion pairing 



Free ions alone are not adequate for treating the low temperature critical region, 
since positive and negative ions will often combine into strongly bound neutral dimers 
or Bjerrum pairsil. This process can be treated as a reversible chemical reaction, say, 
(+) + (— ) # (+,—), leading to equilibrium densities of free ions and dipoles, varying 
with T and In a continuum model, however, there arises a serious question as to 

precisely what configurations are to be considered as bound pairs0'il'0. In practice, this 
relates directly to the problem of determining the proper association constant K(T). Bjer- 
rum's original approachU was to introduce a temperature-dependent cutoff distance that 
would represent, in some sense, the size of dipolar pair. Later Ebeling0, using systematic 
cluster expansions, obtained more elaborate expression for K(T); it turns out, however, that 
Bjerrum's form is reproduced asymptotically to all orders at low temperatures. But for a 
lattice system the situation is intrinsically simpler because a clear and acceptable definition 
of a bound ion pair is two oppositively charged ions occupying neighboring sites. (Pairs 
separated by further distances, second nearest neighbors, etc., may be regarded as distinct 
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Table II: Coexistence curve parameters for 3-D cubic lattices according to pure DH theory; HC 
denotes the continuum hard sphere system, i.e., the RPM 



model 


1 C 


Pi 




sc 


0.101767 


0.007869 


0.0996 


bcc 


0.100617 


0.005908 


0.0759 


fee 


0.096637 


0.004755 


0.0596 


HC 


0.061912 


0.004582 


1 



species and could be considered separately, if necessary^.) 

Following this convention, we introduce the density p 2 and chemical potential fi 2 of Bjer- 
rum pairs which we suppose, initially, behave like ideal lattice particles. The condition of 
chemical equilibrium, fi 2 = 2/ii, ensures that p* and p 2 are interconnected via the law of 
mass action. To this end, let 

Zl = Af/^oCiK 1 , z 2 = Al/(v 2 ( 2 )e^ (30) 

denote activities of free ions and pairs, respectively, where the Aj denote the de Broglie 
wavelengths for which we have A + = A_ = Ai and Ai = A 2 (see Ref.@) while Ci and C2 
represent the corresponding internal configurational partition functions. In terms of the 
activities, the law of mass action states z 2 = \Kz\ 1 from which follows^ 

k ( t ) = = un (si) 



A| VCi 

This definition of K as the internal partition function of a dipolar pair leads naturally to 
the basic expression 

K(T) =v J2 e-^ (a "" ;T) = t; c e-^ (0;T) , (32) 

tin 

where <p(0; T) is given by ([lOD . 

By using the potential-distribution theoremB, we can then write the free ion density as 

pi = Zie -*/ kT = vvCte-^e-^^/kl (33) 

where ^> is the potential of mean-field force and Jxf 1 is given by ([14]) (with d = 3) since neutral 
particles do not contribute to the electrostatic interactions. The second exponential factor 
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here accounts for all the non-Coulombic interactions, since the ionic terms are already taken 
into account by the factor with exponent pif . Hence, only a hard-core factor is required: 
this may be taken as the probability that a given lattice is empty, namely, 1 — p\ — 2p 2 . In 
total the ionic chemical potential may thus be expressed as 



Pi = In 



Pi 



+ In 



A? 



(34) 



\-p\-2p\) \VqCi, 
To obtain a complementary expression for p 2 we appeal to the Bethe approximationS It 
corresponds to the zeroth-order term in the series expansion of the grand-partition function 
for dimers with no attractive interactions and yields 

(2p*/c )[l-(2p*/c )] 



Z2V0 



Thence we obtain 
p 2 = 2 In 



2P2 



Pi - 



+ In (2p*/c ) + ln[l - (2p* 2 /c )] - In 



4A^ 



(35) 



(36) 



\-p\-2p\) v rzi u/ v rzi U/J VC 2 ^o 

On using the law of mass action, the Bethe approximation also yields an equation for p 2 , 
namely, 

(2p*/c )[l-(2p* 2 /c )} = ( p\ \ 2 co &2ll E i_ pqm 

[\-pl-2plY ~ \i- pl -2p* 2 ) 4 ■ 1 ] 

Taking into account that the dimer density should increase as the free-ion density increases, 

we may solve to obtain 



Pi 



Co 



1 - c pl exp 



2vra 3 
3v T* 



P 



6 



x 2 + 6 



1/2" 



(38) 



Since the dimers are neutral, they do not add to the DH interaction energy which retains 
the form (|T3|) . For the total free energy we then have 



f = f M + f 1 
in which we recall that 



2f Id Cm) + f Id (P2) + f El (Pi), 



(39) 



(40) 



x = K,a = Ana pl/v T*. 

Now we may note that the free energy density can be found by integration of (|34]) or (|36] 
with respect to pi or p 2 , respectively. Comparing the resulting expressions yields 



Vl d {\pi)^ = -Pi InpJ -{I -Pi- 2p\) ln(l -pi- 2p\) - pi In (A?/v ) 



(41) 



13 



which can be obtained independently by noting, that the free volume available for an ion is 
proportional to 1 — p\ — 2p\ (see also Ref.H). In addition we get 

f 2 Id v = - 2p* 2 hypl - p* ln(2 P ;/c ) + (c /2) ln[l - (2p* 2 /c )] 

- p\ ln(l - 2 P ;/c ) - p\ In {A[/( 2 v Q ) . (42) 

The equation of state for DH theory with ideal dimers then follows from (p. As in the 
continuum RPM0, one finds that the lattice DHBj theory merely superimposes the pressure 
of an ideal lattice gas of Bjerrum pairs on the DH pressure for the free ions. Additionally, 
the ideal-gas term for the free ions is changed somewhat, since the hard-core interactions 
intrinsic to a lattice gas, appear in the entropic contributions to the free energy via the 
fraction of available lattice sites. However, this does not affect the critical temperature of 
the liquid-gas transition, since it is primarily the free ions density that governs the phase 
separation. 

One now finds, just as in the continuum modelS, that the coexistence curves for all the 
lattices have a banana-like shape: see Fig. 2. This is simply a consequence of the rapid 
growth of the number of neutral dimers as the temperature increases. Indeed, the overall 
critical density is predicted to increase by a factor of 5.02 for the sc, 6.38 for the bcc and 
10.25 for the fee lattice, taking the values p* = 0.03982, 0.03769, 0.04878, respectively. Since 
p* lc does not change when adding ideal dimers, and decreases when the lattice symmetry is 
enlarged, the fact that the overall critical density for fee lattice is greater than for the bcc 
and sc lattices is surprising; but, no doubt, the increased coordination number, c , serves to 
enhance the formation of dimers. 



C. Dipole-ion interactions 

The banana-like shape of the DHBj coexistence curves is clearly unphysicafi Indeed, 
as noted by Fisher and LevinS, the next terms in the expansion of the free energy with 
respect to density take into account the effects of screening of the bare dipole field of a 
Bjerrum dimer by the free ions. As shown in Ref.0, this solvation effect reduces the free 
energy of an ion pair in the electrolyte. It also eliminates the unphysical banana form of the 
coexistence curve and produces better agreement between the critical point predictions and 
the estimates from simulations. 
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Figure 2: Predicted phase diagrams of gas-liquid coexistence: (a) sc lattice with inclusion of Bjer- 
rum pairing alone; for bcc and fee lattices the coexistence curves have a similar form. Coexistence 
curves for (b) sc, (c) bcc, (d) fee lattices based on the full DHBjDI theory. 

Proceeding in this direction, the dipole solvation energy f DI can be calculated via the 
standard DH charging procedure^. Moreover, for the lattice case it turns out that one can go 
substantially farther than for continuum-space models. Indeed, under a few very reasonable 
approximations, we can obtain closed analytical expressions. Consider the positive ion of a 
dipolar pair. Instead of (H) we now have 



where the prime on the summation means that the site with the negative ion is excluded. 
Owing to the symmetry, the potentials of the negative and positive ions differ only in sign, 
the energies being equal. Hence we obtain 



where the potentials (p(sL nn ) can be calculated using the DH expression (^j) separately for 
the contributions arising from the negative and positive ion of each pair. After some algebra 




(43) 




(44) 



nn 
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the results for the cubic lattices can be written in the general form 



V>+ = <p+(x) - ip + (x = 0) 

2nqa 2 1 f f 1 + (c - 1) J(k) - c J 2 (k) 



(45) 



1 



3L>v c + 1 i J k 1 + x 2 /6 - J(k) 
with appropriate values of Cq and «/(k) for each of the lattices: see Eq.(^3|). Utilizing the 
definition of -P(C) i n © enables us to rewrite this in the more convenient form 

„2> 



in which 



G{x 2 



2n( l a2 1 ri r T 2 , rr 



x 2 [cqX 2 + 6(c + 1)] 



P 



6 



6(x 2 + 6) \x 2 + 6 / 

Then the DH charging process may be implemented straightforwardly to yield 

1 



./ 



7di Va 2 1 



3Dv 2 c + 1 



12 x 2 J Q 



(46) 



(47) 



(48) 



with the corresponding chemical potentials 

4tt 2 



ft? 1 



a 6 p* 2 



3(c + l)v 2 T* 



Co 1 
12 x 4 



G{x 2 )d{x 2 ) - G{x 2 ) 



P? 1 



7T 



3(c +l)v 2 T* 



and pressure 



P DI v 



f DJ v + Pi' 1 fi + i% 1 & 



t-.DI * i -DI * 



(49) 



(50) 



(51) 



The only matter not yet taken into account is that, owing to dipole-ion interactions, the 
excess chemical potential will appear also in the law of mass-action. Since at the densities 
of interest for criticality we suppose Bjerrum pairs interact only with free ions — in the 
continuum-space RPM dipole-dipole interactions appear in the next higher term in the 
series expansion — the Bethe approximation for the dimer activity remains adequate. Thus 
we obtain an equation, defining implicitly the pair density function of p*, namely, 

6 



(2p* 2 /co) [1 - (2p*/c )] = ^ c o(Pi) 2 exp{|^ 



P 



x 2 + 6 



(52) 



+2/2?V,p*)-/2^ (pl,pl 



DI/ * 
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Table III: Critical parameters predicted by the full DHBjDI theory. For comparison, values for the 
RPM are also given. 



model 


± C 


Pi 


sc lattice 


0.09666 


0.03041 


bcc lattice 


0.08931 


0.02563 


fee lattice 


0.08064 


0.02708 


RPM: DHBjDI 


0.0554-0.0522 


0.0244-0.0259 


RPM: simulations 


0.049 


0.06-0.085 



This completes the principal task and allows the construction of coexistence curves: these 
are shown in Fig. 2. Clearly, in the lattice models the dipole-ion interactions are also 
crucial to repair the unphysical banana-like form produced by Bjerrum association alone. 
The numerical estimates are presented in Table III. For comparison, the predictions of 
continuum-space DHBjDI theory and of the RPM simulations are also listed. The predicted 
lattice critical temperatures are now 1.5-1.65 times greater than the value given by the 
simulations and the theoretical results of Levin and FisherM The critical densities, however, 
are quite close to the continuum model predictions, but are significantly lower than the 
simulationsill§lll. 

IV. SUBLATTICE ORDERING 

So far we have dealt only with an intrinsically low-density picture of the system. Our 
description of the dense phases, although partially represented by the liquid side of the co- 
existence curve, has been seriously incomplete. On the other hand, lattice theories provide 
a particularly natural first approach to studying ordering in solid phases. Indeed, the ques- 
tion of principal interest for us will be the possibility of ordering similar to that observed in 
ionic crystals. We will, in fact, find that a DH-based theory yields a phase diagram with no 
gas-liquid criticality but, rather a tricritical pointffli. 
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A. DH Mean Field Theory 

Let us start by considering a general rf-dimensional bipartite lattice, that can be divided 
into two sublattices of the same form. Suppose, Nj[, iV7_ and jV#, Ng are the numbers 
of positive and negative ions on sublattice A and B, respectively, subject to the neutrality 
constraint Nj[ + N^ = Ng + Ng = N/2. Consider the sublattice with an excess of positive 
ions (say, "A" for the defmiteness) , and define the corresponding order parameter by 

y= N j- Nl - (53) 
This will have a vanishing mean value in a disordered phase but will be positive in an 
appropriately ordered phase. 

The entropic part of the free energy density corresponds to ideal ions and is thus now 
given by 

/" = -p\ lnp* - (1 - p*) ln(l - pi) - \p\ [(1 + y) ln(l + y) + (1 - y) ln(l - y) - 2 In 2] . 

(54) 

To estimate the electrostatic part of the free energy, the extended DH approacb0 suggests 
that we begin with an inhomogeneous version of Poisson's equation for the potential at a 
general site r due to all the ions when an ion of type a is fixed at the origin: this states 

DAip a (r) = -C d J2 qrpr(r)g T A r ) + C d q a 5(v), (55) 

r 

where p T (r) = p*(r)/vo is the bulk density of ions of species r while g T ,a{ r ) is the ion- ion 
correlation function. Approximating the correlation functions by simple Boltzmann factors 
and then linearizing provides a DH equation. 

However, we must now allow for an overall nonzero charge density on each sublattice 
given by 

UA = N/2 = PlVq ' UB = N/2 = ~ PlV<1 ' ^ ' 

These charge densities generate an additional "background" potential, $(r), which does not 
contribute to the correlation functions since it is independent of what type of charge is 
placed at the origin. For sublattice A we thus have a linearized Poisson-Boltzmann or DH 
equation 

= -^p* m = (57) 
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Recalling the definition of the lattice Laplacian (|5]), and taking into account the symmetry 
between the sublattices we have 

$(r A ) = -$(r B = r A + a nn ), (58) 

and conclude that the background potentials are 

*<"> = = m- < 59) 

If we now put (f>A = ¥>a — &(ya) and accept the approximation g T ,o-(r) — exp(— (3(p) in 
(^5|), the equation for the local induced potential if a reduces to the lattice version of pure 
DH theory (|). This reflects the electrostatic superposition principle, i.e., the total potential 
is simply the sum of a "background" potential due to nonzero average charge density and 
the DH screening potential so that 

CdQ a 2 f e ik r , , 

+ + (60) 

Now, following the DH approach combined with a mean-field description of the ordering, we 
find the potential if) due to all ions except the one fixed at the origin to be 

<p A {r A = 0) = - Y] [$(a nn ) + v DH (*nn)] , (61) 
ipA = MO) -<Pa(0)\x=o, (62) 
which, on taking into acount (|5B]) and (|5Tf), yields 

^ = -A + *°' (63) 



with if> DH given by the same expression as if>i in (p!2|). 

Finally, the DH charging process gives the total free energy density (of both sublattices) 
as /= f° rd + f DH where f DH = f Id + f El follows from (@) and (|TJ), while 



f° rd = g^g {P1 1 Y - ^ [(1 + V) ln(l + V) + (1 - y) ln(l - y) - 21n2] . (64) 



Note that this result implies that the electrostatic part of the ordering energy is negative 
(/ = —F/kTV) as it should be since it describes the interactions between charges of opposite 
signs. The ordering term also yields additions to the chemical potential and pressure, namely, 

fi° rd = + + V) ln(1 + V) + (1 " V) ln(1 " V) ~ 2 ln 2] ' (65) 
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frd = Cd a d pfy 2 /8dv 2 T*. (66) 

Now the possibility of sublattice ordering is explored by seeking minima of f° rd with 
y 7^ 0. This leads to 

= h, ( ±±*) . (67) 
Expanding for small y in the standard way yields the solution 

Z/-V3(|-l) 1/2 , (68) 

in which the A-line, along which second order phase transitions occur, is given by 

p* x (T) = (Adv /C d a d )T*. (69) 

The simplest way to find the anticipated tricritical point is to consider the intersection of 
the spinodal with the A-line^. This can be found by computing dp(p* 7 y) / dpi with y defined 
by (^), equating to zero and setting p* = p\(T) after taking the derivative. A more readily 
justifiable, but also somewhat more sophisticated procedure is to study the stability matrix 
for the free energy, which is now a function of two order parameters, namely, y and p\ . Both 
methods lead to the same equation for tricritical point, which reads 

AddPf 2d 
p* tr dx 2 \x 2 + 2d 



+ — = 0, (70) 

1 — n* 9 

Ad 1 Ptr z 



where P(C) is the lattice Green's function 



B. Results and Discussion 



In d = 3 dimensions simple cubic and body centered cubic lattices are bipartite and 
sublattice ordering is possible. Note, that the calculations presented above did not use any 
extra properties of lattice symmetry. Indeed, the electrostatic interaction energy of the 
"charged" sublattices depends only on the excess charge density, that is on y, and is thus 
the same for both lattices. This is also true as regards the entropy of sublattice ordering. 
We find that the tricritical point parameters are 

T t * ri = 0.3822 (sc), 0.4865 (bcc), (71) 
p\ H = 0.3649 (sc), 0.3576 (bcc), (72) 




Figure 3: Phase diagrams for ionic lattice systems with sublattice ordering: (a) sc, (b) bcc. Also 
shown by broken lines are the gas-liquid separation curves predicted by pure DH theory for (c) the 
sc and (d) the bcc lattice. 



while the A-lines may be written 

T A */p* = |tt ~ 1.047 (sc), ~ 1.360 (bcc). (73) 

The full predicted phase diagrams are presented in Fig. 3. 

It is instructive to compare our results with previous simulations for the sc latticeS@. 
These yield T 4 * ~ 0.14 — 0.15 which is only about 40% of our theoretical estimates (|7T|). On 
the other hand, for the tricritical density the higher estimate p* tri ~ 0.48 of Ref.0 is probably 
more reliable than p* tri ~ 0.38 of Ref.@ (which compares rather well with our theoretical 
values), since the former simulations used larger lattice sizes and computed more points on 
the coexistence curve. It must also be noted, however, that both these simulations employed 
the discretized continuum or 1/r Coulombic potential in place of the lattice form we have 
used. 

Fitting a straight line at low T to the A-line data of Stell and Dickmanil yields a slope 
Tjjy ' p\ ~ 0.6 which may be compared with our value of 1.047. (For another comparison 
one might note that the generalized DH theory for the continuum RPM0 predicts damped 
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charge- density oscillations setting in on a locus T^/p* K ~ 0.3 while undamped oscillations 
are predicted beyond the locus T^/p* x ~ 9.) 

Also of interest is the value of T£(p*) at close packing, i.e., p* = 1. For the sc lattice 
our treatment predicts T x (l) ~ 1.047 which, by virtue of the mean-field character of the 
theory, is likely to be a significant overestimate. Indeed, extrapolation of the Dickman-Stell 
data suggests 7^(1) ~ 0.6 while Almarza and Enciso0 obtain T^(l) ~ 0.515; but, again, 
these simulations employ a discretized 1/r potential. Our analysis indicates higher values 
of T t *j and T^(l) for the bcc lattice than for the sc. But, perhaps, surprisingly, this is just 
the opposite of what Almarza and Enciso find§. In addition to using the lattice Coulomb 
potential, our treatment at this point has neglected the formation of Bjerrum ion pairs. 
In fact, it seems quite feasible to include pairing in the theory along with allowance for 
sublattice ordering (since the ions of a dimer pair will reside on complementary sublattices). 
It is possible that this improvement of the theory will result in lower transition temperatures 
for the bcc lattice relative to the sc lattice. 

Previous theoretical discussions of the sc lattice RPM have been presented by Stell and 
colleagues0'0. In an initial mean field approach^, the long-range Coulomb potential (taken 
in discretized 1/r form) was first reduced to an effective nearest-neighbor interaction. The 
value of the tricritical temperature, T* ri = 2, obtained by direct mean-field lattice calcu- 
lations, was then scaled down by a factor derived by comparing energy magnitudes. This 
approach suggested T* ri ~ 0.3 and p* tri = |, the latter value being merely a consequence of 
using a nearest-neighbor mean-field approximation. More recently, Ciach and Stelli have 



adopted a single-ion lattice potential, as, in fact, given by (||) and (jlQf) . This corresponds 
more closely to our treatment but they entirely neglect the cooperative screening which must 
occur and which is included in our DH-based treatment. (Note that at higher densities the 
screening effectively takes place via "holes" in the ordered or close-to ordered lattice charge 
configurations). The new treatment^ reproduces p* tri = ~ (for the previous reasons) but 
gives T* ri ~ 0.6, which is worse than the previous result as compared with the simulations; 
however, no energy rescaling is now performed. 

As we have seen, both by our own theoretical analysis and through the simulations, 
the sc and bcc pure Coulomb lattice systems display no gas-liquid phase separation as 
such. Indeed, in Fig. 3 we have plotted the coexistence curves for the two lattices that 
are predicted by pure lattice DH theory with no allowance for the possibility of sublattice 
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ordering. Evidently, these coexistence curves lie entirely within the two-phase coexistence 
region of dilute disordered vapor and the high density ordered "crystal". Consequently, 
the order- disorder transition entirely suppresses gas-liquid criticality in the simplest discrete 
ionic system with only single-site hard cores. If, instead, the hard-core repulsions extend over 
more lattice sites — or, equivalently, if a finer level of spatial discretization is employed (as, 
of course, is more-realistic for continuum systems) — then, as revealed by simulations0'0, 
a normal gas-liquid transition and critical point is restored. At the same time, ordered, 
crystal-like phases appear only at relatively higher densities as characteristic of real solids. 

While a DH (or, even, a DHBjDI) theory might be attempted for a more finely discretized 
model, the clearly evident complications do not make this a promising prospect. On the 
other hand, by adding to a purely ionic lattice system strong short-range attractive potentials 
(say, designed to represent neutral solvent properties^), more elaborate phase diagrams can 
be anticipated. Indeed, by approximations that again neglect all screening effects, systems 
displaying both a tricritical point and a normal critical point have been obtained^. Our 
more complete treatment could readily be extended in the same spirit. 



V. CONCLUSIONS 



By solving exactly lattice versions of the usual Debye-Hiickel equation, we have derived 
closed expressions for the free energy of general, <i-dimensional ionic lattice systems with 
single-site hard core repulsions. For d > 3, gas-liquid transitions are predicted at low 
temperatures and densities. As in the corresponding DH-based theory for the continuum 
restricted primitive mo improvement of the theory at low temperatures demands both 
allowance for (+, — ) ion pairing, to form nearest-neighbor dipolar dimers, and the solvation 
of the resulting dipolar pairs by the residual free ions. The predicted critical temperatures 
for the sc, bcc and fee lattices in d = 3 dimensions then lie 60 — 70% higher than given by 
continuum DH-based theories; but the critical densities are relatively closer. These results 
accord with the general tendency of lattice theories to overestimate the stability of the 
corresponding low-temperature continuum phases. 

At higher densities in a lattice theory it is imperative to allow for sublattice ordering of 
the positive and negative ions. By introducing an appropriate order parameter we have ex- 
tended analysis to treat general, (^-dimensional bipartite lattices at a combined Debye-Hiickel 
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and mean-field ordering level. Our unified theory yields, in an accord with recent lattice 
simulations, a complete suppressiion of gas-liquid phase separation and criticality by order- 
disorder transitions that occur at higher temperatures. At high densities and temperatures 
a classical second-order A-line is predicted; but this terminates at a tricritical point at a den- 
sity, for the sc and bcc lattices, p\ ri j ' p* max — 0.36 and a temperature, Ttril^max — 0-4 ~ 0.5. 
At lower temperatures the first-order transition is from an exponentially dilute vapor to an 
almost close-packed ordered ionic lattice. 

Our treatment can be extended in various directions. Indeed, there are preliminary 
indications that by considering strongly anisotropic three-dimensional lattices, gas-liquid 
separation may be restored, possibly, together with distinct order- disorder transitions. It 
is relevant to note in this connection that DH theory for continuum ionic systems predicts 
increasing values of gas-liquid critical temperatures when the dimensionality is decreased^!. 
Thus lattice anisotropy might mimic lower dimensionality. 

Although the direct applicability of our results to ionic systems is clearly limited, we feel 
the approach developed in Sec. IV may prove helpful in describing the behavior of defects 



potential are desirable and might cast some light on the role of short-range interactions and 
geometric constraints in strongly coupled ionic systems. 
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